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Large Eddy Simulations (LES) of gyrokinetic plasma turbulence are investigated as interesting candidates 
to decrease the computational cost. A dynamic procedure is implemented in the GENE code, allowing 
for dynamic optimization of the free parameters of the LES models (setting the amplitudes of dissipative 
terms). Employing such LES methods, one recovers the free energy and heat flux spectra obtained from 
highly resolved Direct Numerical Simulations (DNS). Systematic comparisons are performed for different 
values of the temperature gradient and magnetic shear, parameters which are of prime importance in Ion 
Temperature Gradient (ITG) driven turbulence. Moreover, the degree of anisotropy of the problem, that 
can vary with parameters, can be adapted dynamically by the method that shows Gyrokinetic Large Eddy 
Simulation (GyroLES) to be a serious candidate to reduce numerical cost of gyrokinetic solvers. 



I. MOTIVATION AND CONTEXT 

In the area of fluid turbulence, theories are usually 
based on the notion of an inertial range in which the en- 
ergy cascades from larger scales to (somewhat) smaller 
scales mediated by the quadratic nonlinearity. The role 
of the smallest scales is then to dissipate energy in the 
so-called dissipative range. In numerical simulations, this 
picture has led to the development of Large Eddy Simu- 
lation (LES) techniques that are based on the idea that 
neglecting the small scales can be compensated by intro- 
ducing a dissipative model for the eddy viscosity^. 

A Direct Numerical Simulation (DNS) is supposed to 
retain all the scales from the injection range down to the 
dissipative range. This requires an enormous numerical 
effort in the case of high Reynolds number flows. On the 
contrary, a LES coarsens the simulation grid and only re- 
tains the largest scales (which are problem-dependent), 
while the small scales (which are assumed to be univer- 
sal) are replaced by a model. In Fourier space, such a 
coarsening can be seen as the action of a low-pass filter. 
Since the scale range is truncated, the dissipation scales 
can not be reached, and the modeling basically consists 
of the introduction of artificial dissipation mechanisms. 
From a more mathematical viewpoint, one notes that the 
filtering operation docs not commute with the nonlin- 
ear term that transfers energy from largest to smallest 
scales, and the major problem of LES consists in finding 
a satisfying closure for representing the influence of the 
unresolved scales. 

Recent gyrokinetic studies have shown that Ion Tem- 
perature Gradient (ITG) driven turbulence exhibits a di- 
rect and local cascade of a nonlinear invariant, namely 
the free energy— Such a cascade is analogous to the 
kinetic energy cascade in three dimensional Navier- 
Stokes turbulence. The important difference is that the 
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quadratic conserved quantity in fluid dynamics is the ki- 
netic energy, while it is the free energy in gyrokinetics. 
The latter quantity is the sum of both the perturbed 
entropy and the electrostatic energy. Transfers between 
entropy and electrostatic energy are ensured by the mag- 
netic curvature and parallel dynamics terms^ 

Thus, adapting LES methods to gyrokinetics is quite 
a natural idea. In a pioneering study, Smith and Ham- 
mett have applied LES techniques to a set of gyrofluid 
equations^ The use of a hyper- viscosity model was found 
to provide better agreement than a simple Smagorinsky- 
type eddy viscosity^ Promising recent comparisons be- 
tween gyrokinetic LES and highly resolved DNS have 
motivated the present work£ As is well known, LES ap- 
proaches have to face two distinct difficulties. First, a 
suitable model has to be designed to mimic the dissi- 
pative effect of the small scales. Second, the free pa- 
rameters of the model have to be determined in order 
to ensure that it creates the correct amount of dissipa- 
tion. In a previous study^, the first difficulty has been 
addressed: the feasibility of Gyrokinetic Large Eddy Sim- 
ulations (GyroLES) has been demonstrated with perpen- 
dicular hyper-diffusion models. However, the second dif- 
ficulty remains to be tackled. In practice, up to now, the 
free parameter setting the amplitude of the dissipative 
term modelling the influence of the neglected scales had 
to be determined via a trial and error process. The main 
objective of the present study is to overcome this prob- 
lem by adapting the dynamic procedure^ to GyroLES. 
The dynamic procedure is an optimization approach that 
allows to calibrate the model amplitude in the course of 
the LES. 

The remainder of the present paper is organized as 
follows. After a brief review of the GyroLES formalism in 
Section HIl the effect of truncating small scales is studied 
in detail, and the dynamic procedure for gyrokinetics is 
discussed in Section Hill Numerical results obtained for 
various logarithmic temperature gradient and magnetic 
shear values are presented in Section HVl 
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II. LES FORMALISM IN GYROKINETICS 

In the following, the nonlinear gyrokinetic equations 
are solved by means of the GENE code^. Although 
a more comprehensive code version including nonlocal 
effects is at hand^, for simplicity we restrict ourselves 
here to the local code version. Only electrostatic fluc- 
tuations are considered, with a fixed background mag- 
netic field Bq and adiabatic electrons. Field aligned 
coordinates are usecU^, with the assumption of circu- 
lar concentric flux surfaces^. The GENE code uses a 
delta-f splitting of the unknown distribution function: 
Fi = Foi + jki with the normalized equilibrium distri- 
bution function F Ql = e - '"!! _AlS ° , where n = miv\/{2B Q ) 
is the ion magnetic moment (mass to,), v± and Uii arc 
respectively the velocity coordinates perpendicular and 
parallel to the magnetic field. Unknowns are Fourier 
transformed along coordinates perpendicular to the mag- 
netic field (x,y) —> (k x ,k y ). The gyrokinetic Vlasov 
equation for ions guiding center distribution function 
fki {kx , ky , 2, v\\ , fi, t) then reads: 



dtfki = L[f kl ] + N[^ k ,f kt ] - D[f kl ], 



(1) 



where L represents linear terms, N the quadratic nonlin- 
carity, and D the numerical dissipation terms. 

The linear terms can be written as L = Lb + Lq + Lii , 
where Lo[f k i] is the drive due to logarithmic density 
and temperature gradients (u) n i and ujti), isj/ti] cor- 
responds to both the curvature and the gradient of the 
magnetic field B (referred to as "curvature" in the fol- 
lowing), and L||[/fej] is the term describing the parallel 
dynamics: 
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Here, hki = fki + Z l F oi J ok (j) k T e o/T l o is the nonadiabatic 
part of the distribution function, with the ions charge 
number Zi and the ion thermal velocity VTi- Tio and 
T e o are, respectively, the ion and electron equilibrium 
temperature, Jo k is the zeroth order Bcssel function cor- 
responding to Fourier transformed gyroaverage operator, 
and <j>k is the electrostatic potential. The two terms K x 
and K y are due to magnetic field curvature and gradient 
introduced by the magnetic geometry^. 

N is the nonlinear term describing the perpendicular 
advection of the distribution function by the E x B drift 
velocity: 

N[(pk,fki] = - y^XKky - k xk' y ) Jok'0k'/(k-k')i > ( 5 ) 



which has the fundamental role of coupling different per- 
pendicular k x and k y modes. 

Numerical dissipation terms in GENE have the general 
form: 



D[f H ] = a x k^f ki + a v k n y f kl + a z dtf kl + a n d 4 (| f kl , (6) 

where the coefficients a x and a y are usually set to zero, 
while a z = 0.1 and a v ^ = 1 have been shown to be well 
adapted in a wide range of cases^. 

The electrostatic potential <j) k is given by the quasi 
neutrality equation: 



i fc > H +^2£[i-r (6, 



= 7tB J dv\\dixJo k fk , 
(7) 

where (4> k ) FS = (J Jdz(f>k) / (J Jdz), stands for the flux 
surface average of the electrostatic potential, To (6,) is 
the modified Bessel function applied to the argument 
hi = VTik\/VL 2 ci . Electrons are assumed adiabatic: n e = 
q e n e0 (4> k - ((f>k) FS ) /T e0 . Since a single gyrokinetic ion 
species is considered, the species indices are omitted in 
the following for the ions distribution function: f k = f k i. 



A. Filtered gyrokinetics 

In a gyrokinetic LES, the most suitable coordinate 
subspace for coarsening the grid is the perpendicular 
wavenumber plane (k x , k y ) since it generally requires 
fairly high resolution. Obviously, the objective of the 
LES technique is to reduce the number of grid points in 
(k Xl ky) space. The coarsening procedure can be imple- 
mented by applying a Fourier low-pass filter, with the 
characteristic length A. The employed cut-off filtering 
has the effect of setting to zero the smallest scales charac- 
terized by all modes larger than f/A, as shown in Fig. m 
If one denotes the action of the filter on the unknowns 
by 7TT , the filtered gyrokinetic equation reads: 



dtfk - L[f k ] + Nfc, f k ] + r SjADNS - D[f k ] , (8) 

where a new term appears from the filtering of the non- 
linear term: 



A,A E 



N[^ k , f k ] - N[^ k J k ] 



(9) 



At this point, it is important to note that Eq. (|9]) is 
the only term which contains the influence of the scales 
A DNS which we want to filter out from (4> k , f k ). Wwe 
will refer to it as sub-grid term in the following. The 
GyroLES then consists of finding a good model replacing 
this term which only depends on the resolved unknowns 
(<j> k , fk), on the characteristic length of the filter A, and 
on some free parameters {c n }. 
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FIG. 1. Schematic view of a Large Eddy Simulation: The 
smallest scales (grey area between dashed-dotted and dotted 
lines) are retained only in a DNS, while they are modeled in 
a LES model; LES only retain the area inside the dashed- 
dotted line; alternatively or additionally, a test filter can be 
used (hatched area, solid line). 

B. Free energy and sub-grid term 

As has been shown both thcoretically-i^r— and 
numericall y 16 ! 17 , the free energy is a relevant quantity 
for studying gyrokinctic turbulence. The free energy is 
defined as: 

= rnOrn yy f ndzd d h-nfu 

with the volume V = ^ *£ ky I dz/B . 

The dynamics of the quantity £ can be derived from 
Eq. ([1]) by the action of the "free energy operator" S on 
the distribution function j k ^\ £ = ^E[f ki ] with 

B[&] = 9ffEE/ -dzdv^^ k . (11) 

One thus obtains: 



is injected at large scales by the background gradients 
and dissipated at various smaller scales by the dissipation 
terms D. It is important to note in this context that the 
parallel advection term (iy), the magnetic term {Lb ), 
and the nonlinear term (TV) have a null contribution to 
the total free energy balance. 

III. DEVELOPING A GYROKINETIC LES MODEL 

As is well known, a naive truncation of small scales 
can lead to a pile-up of free energy at the smallest 
scales which are retained in the filtered simulation^ A 
good LES model is thus required to dissipate the correct 
amount of free energy. In the following, the role of sub- 
grid terms in the free energy balance will be studied in 
detail. A model will then be developed which agrees as 
much as possible with the desired sub-grid properties. 

A. Sub-grid term and dissipation of free energy 

The nonlinear term has the fundamental role of trans- 
ferring free energy across perpendicular scales without 
affecting the global free energy balance. This property is 
expressed by 

E[N[<j> k J k ]]=Q, (14) 

simply reflecting the fact that the nonlinearity has a Pois- 
son bracket structure and, consequently, its integration 
cancels upon integration. For the same reason, if a filter 
is introduced, the following property holds: 

E[N[fa,Tk}]=0, (15) 

where E is the filtered free energy operator defined in 
the filtered space. On the contrary, the filtered free en- 
ergy operator has a non vanishing contribution when it 
is applied to the sub-grid term: 

7^ A dns = H[T SjA dns] = 3 \N[<j>k,fk] ~ N[(j> k , f k }] 

= E{N[<f> k ,f k }} . (16) 

The filtered free energy balance can then be expressed 

as 



d t £ = g-V, (12) 

with the definitions 

G = E [L G [f ki \] , V = E[D[f ki }}. (13) 

This balance is of particular relevance for the design 
of a good model. As pointed out in Rcf.— , Eq. (fT2|) 
involves only quantities which are quadratic in the dis- 
tribution function, like the kinetic energy in fluid turbu- 
lence. Moreover, like the latter quantity, the free energy 



d t £ = Q + 7^ :ADN s - V , (17) 

where filtered quantities are obtained from the action of 
the filtered free energy operator E on the filtered gyroki- 
netic equation (JU). 

Recalling that the free energy is assumed to be injected 
at large scales, then transferred to smaller scales and dis- 
sipated there, one can expect that the sub-grid contribu- 
tion to free energy balance (fT6| will be negative. Indeed, 
one can expect that a large majority of the free energy 
injection will not be affected by the filtering: Q « Q. It 
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FIG. 2. Contribution of the sub-grid term to the free energy 
balance as a function of time, for different test-filter widths 
A. 



follows that the DNS dissipation can be approximated by 

T) r=i T) 7^ ^DNS ■ 

The time evolution of the sub-grid contribution to the 
filtered free energy balance (CBC for ITG range of pa- 
rameters) for different values of the filter width A is 
shown in Fig. [2j One can see that the sub-grid con- 
tribution is always negative, implying that the sub-grid 
scales act as a free energy sink, like it is supposed to£ 
More precisely, one observes that the amplitude of the 
dissipation ensured by the sub-grid scales increases with 
the filter width. This means that a model should behave 
like 



the inertial range. Based on the recent finding that ITG 
turbulence also exhibits a local and direct cascade of free 
energy^, we assume, in close analogy, that the free en- 
ergy flux S£ is constant from scale to scale in the (k x , 
k y ) plane perpendicular to the magnetic field. The free 
energy has the dimension of an energy density, so that 
the free energy flux eg is an energy density per time, 



r 



where r and £ represent characteristic time and length 
scales. It is reasonable to assume that the model depends 
only on the free energy flux e$ and the filter width A, 



M 



A a e F k n h k . 



Moreover, from dimensional analysis we know that 
[M] = T- l [h k ], so that [3 = 1/3 and a = n + 1/3. The 
last relation allows to fix the unknown filter width expo- 
nent a accordingly to the model parameter n. The model 
thus becomes 



M 



'c.A: +1/3 fc: 



■ Cyl\ y Ky 



I'M 



(20) 



Since the derivative order n is positive, the filter width 
exponent a = n + 1/3 is also positive, in line with the 
numerical results in the previous section. Moreover, the 
model coefficients are dimensionally related to the con- 
stant free energy flux across scales via [c x ] = [c y ] — 



1/3 



It is interesting to note here that the model co- 



efficients are constants, just like the free energy flux. 



M(c,AJ k ) = A a M'(c,f k ) 



L A ADNS 



C. Dynamic procedure for gyrokinetics 



B. A model for sub-grid scales 

A simple dissipative model for GyroLES which has al- 
ready been used previously^ is given by 

M(c ± J kl ) = c ± kih kl . (18) 

The optimal value of c± for the CBC parameters can be 
found, e.g., through trial and error. However, this model 
is not taking into account the filter width dependency A 
observed in the previous section. Moreover, the use of 
k± implies that the relative dissipation in k x and k y is 
fixed. A more flexible model which takes into account the 
anisotropy (c x and c y ) and the filter width dependency 
(A^y) is given by 

M=(Xc x K + A^Cyk;)h ki . (19) 

In fluid turbulence, it is common to assume that the 
kinetic energy flux from scale to scale is a constant in 



The dynamic procedure is based on the introduction 
of an additional filter denoted by — and referred to as 
the test-filter. It is characterized by a filter width A that 
corresponds to a "very coarse" grid: A > A > A DNS . 
The gyrokinetic equation associated to the test-filter grid 
can be obtained by test-filtering the gyrokinetic equation 
expressed in the DNS domain: 

dtfk = L[f k ] + Nfa, f k ] - D[f k ] + T K ADNS . (21) 

This equation is equivalent to the LES filtered Eq. © 
with the LES width (A) replaced by the test-filter one 

(A). 

Alternatively, the equation in the test-filter domain can 
be obtained by test-filtering (A) the gyrokinetic equation 
expressed in the LES domain, Eq. (|5J|, 

dtfk = L[f k ] + N[J~ k , A] - D[f k ] + f SjADNS , (22) 
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where we have used the very important property 7TT = — 
of Fourier cutoff filters. Comparing Eqs. ([2Tjl and ([22]) . 
one obtains the Germano identity, 

7S jA DNS = % A DNS + N[4>k,fk] - N[4> k ,f k ] , 

— ^a,a dns + ^k,A ■ (23) 

During an LES, the sub- grid term ^ can be com- 
puted exactly, since it involves test filtering (A) of the 
LES-rcsolved quantities (A). On the other hand, the two 
other terms involve the non- resolved DNS scales (A DNS ) 
and therefore have to be approximated by the model: 

^A,A DNS ~ i %,A DNS (24) 

The dynamic procedure consists of introducing the 
model approximations, Eq. (|24|). into the Germano iden- 
tity, Eq. (f23]> . to obtain 



Thus, these two free parameters of the model can 
be computed dynamically during a numerical simulation 
from Eqs. (|28j) and (|29|) . Since they are always positive, 
it is guaranteed that the model has a dissipative effect 
on the free energy. 

IV. NUMERICAL RESULTS 

In the following, we will present numerical results ob- 
tained by means of the dynamic procedure with the 
GENE code. The set of parameters corresponds to the 
Cyclone Base Case commonly used for studying Ion Tem- 
perature Gradient (ITG) driven turbulence^. Consid- 
ering a minor radius ro/Ro = 0.18, the density and 
temperature gradients are, respectively, tu n i = 2.22 and 
= 6.96, the magnetic configuration is characterized 
by the safety factor q = 1.4 and the magnetic shear 
s = 0.796, with ions and electrons such that T e o/Xio = 1 
and Zi = 1. 



M^M-z + T^. (25) 

Since the model is an approximation of the sub-grid 
term, Eq. Q23|) can only be approximated during an LES. 
Now, one can define the squared distance d 2 which is to 
minimize via 

^ = ((^ a ,a + ^a-^a) 2 ) a . (26) 

where (• ■ • ) A stand for integration over the entire phase 
space. 

As was shown in Sec. MI B[ the model coefficients c x 
and c y can be assumed to be constant in the gyrokinetic 
"inertial range." So provided that the range between 
test-filter and LES scales belongs to this "inertial range," 
the coefficients do not depend on the filter widths (A, A). 

Using Eq. (j2TJ)) , the squared distance can be expressed 
in terms of the model amplitudes c x and c y according to 

d 2 = ^ (t-.a + c * m x + c v m v) ^ > (27) 

where the notations m XtV = ^A" y — A" ^ k" y hk have 

been introduced. 

An optimization of this difference with respect to the 
unknown parameters (dd 2 /dc x = and dd 2 /dc y = 0) 
leads to the expressions 

(m»T s 3 ) (ml) - (myT^) {m y m x ) A 

c x = - — g " T^T (28) 

(m x m y ) A - (to2) a (mg) A 

( m y T A,A/ A ( m l) A - ( TO ^ T A,A/ A ( m v m ^A 

Cy = ± ^ 2 i (29) 

(m x m y ) A - K) A (m 2 y ) A 



A. Nonlinear Gyrokinetic Large Eddy Simulation: Cyclone 
Base Case 

For the reference DNS, a perpendicular grid of N x x 
N y = 128 x 64 is used. This grid has been used both with 
and without a LES model, and the results obtained have 
not been affected, indicating that the simulation is well 
resolved. On the other hand, a minimal perpendicular 
grid for GyroLES should be N x xN v =48x 32, allowing 
the dynamic procedure to work. Indeed, the use of the 
latter involves the introduction (in the LES domain A) 
of a test filter corresponding to a coarser grid, A > A. 
However, it is necessary for the dynamic procedure that 
the domain of the LES grid which is neglected by the test 
filter belongs to the gyrokinetic "inertial" range, so that 
the model coefficients have the same values in the two 
domains. Here, we will employ a test filter width which 
corresponds to the half of the LES domain: A T ~ 2A X , 
A y = 2A X . This means that the optimization in the 
dynamic procedure is related to a sub-grid term 
defined by 24 k x modes and 16 k y modes. The parameters 
given in Eqs. (f2"8"]) and (|2"5|) are computed at each time 
step of the simulation. The parallel and velocity grids 
are kept fixed at N z = 32, N n = 64, and = 8. The 
model order is chosen to be n — 4, leading to 

M4 = (r, A'" + c y A y 3/3 kl) h kl . (30) 

We note in passing that the case n = 2 has also been 
tested; in that case, the obtained results showed a high 
volatility, though. A similar result has been obtained 
by Smith and Hammett^, where hyper-viscosity models 
were found to perform better than viscosity models for 
gyrofluid turbulence. For comparison, we also show sim- 
ulations without any model (LES MO), with a perpendic- 
ular grid identical to the LES M4 one. 
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FIG. 3. Free energy injection spectra (Q kx at top, Q ky at bot- 
tom) for the fourth-order model (M4) at reduced resolution, 
compared with a highly resolved DNS and the case without a 
model (M0). 



FIG. 4. Free energy spectra (£ fca= at top, £ ky at bottom) for 
the fourth-order model (M4) at reduced resolution, compared 
with a highly resolved DNS and the case without a model 
(M0). 



Since it is proportional to the total heat flux Q, the free 
energy injection term Q = LJTiQ is of special relevance for 
comparisons with experimental results and earlier works. 
The comparisons are based on two dimensional wavenum- 
ber spectra of the free energy injection rate, 

= W I ■ < 3i » 

It is understood that this quantity is averaged during 
the quasistationary turbulent state over sufficiently long 
time windows (at least 2000 Ro/vri)- The reduction to 
a one-dimensional spectrum is then simply provided by 

gk x _ gk^.ky ^ gk y _ V"^ gk^.k y ^ 

Comparisons of free energy injection spectra from 
DNS, GyroLES, and a simulation without a model are 



shown in Fig. [3] The k x spectra for all three cases are 
found to exhibit a surprisingly good level of agreement, 
but the k y spectra illustrate that the use of a LES model 
diminishes the accumulation at the smallest scales, im- 
proving the agreement at the largest scales with the ref- 
erence DNS. 

One-dimensional free energy spectra, 8 fcx and £ ky , can 
be constructed in analogy with Eq. (p?2"j) . As can be ob- 
served in Fig. 21 the GyroLES clearly prevents the accu- 
mulation of free energy at the smallest scales. However, 
there still exists an overestimation of the free energy at 
the largest scales when compared with the reference high- 
resolution DNS. 

Since the LES spectra are truncated with respect to the 
DNS ones, estimates for the neglected parts of the spectra 
are required for computing total (integral) values of the 
heat flux and the free energy. Such estimates may be 
provided via a power law regression of the spectral The 
estimate from the GyroLES run yields £ Mi = 1.49 £ DNS , 
while for the case without a model (M0), one obtains 
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£ M0 = 2.36£ DNS . The total heat flux levels from DNS 
and LES are in very good agreement, <2 M4 = 1.06 Q DNS . 
In the case without a model, one finds Q M " = 1.20 Q DNS 
due to an overestimate at the smallest scales. 

For the sake of clarity, all further comparisons will fo- 
cus on the spectra £ ky and Q ky which have been found to 
be most sensitive. The LES model uses n = 4 and filter 
widths such that A = 2A; the perpendicular grid size is 
N x xN y =48x 32. 



B. Robustness while varying the temperature gradient 

As is well known, the logarithmic temperature gradi- 
ent LOTi is a key parameter for ITG turbulence, given 
that the equilibrium temperature profile acts as a source 
of free energy for the system. In the following, the ro- 
bustness of the LES approach is tested for two values of 
the temperature gradient which differ from the nominal 
value; these correspond to a weakly driven turbulence 
case (ujtz = 6.0) and to a strongly driven turbulence case 
(u Tl = 8.0). 

The case of weakly driven ITG turbulence is shown in 
Fig. El The M4 model yields a very reasonable agreement 
with the DNS regarding both the free energy spectrum 
£ ky and the free energy injection spectrum Q ky . The 



total values £ " 



1.02 £ DNS and Q K 



1.25 Q DNS are 



also in good agreement. Without a model, one obtains 
g M = 1,79 £dns and quo = i 04 Q DNS . The latter result 

is accidental, however, and results from a compensation 
between an underestimation at large scales and an over- 
estimation at small ones. 

Fig.|6]displays the results for the case of strongly driven 
ITG turbulence. The LES is found to systematically 
overestimate the DNS free energy spectrum £ ky , while 
the prediction of the free energy injection spectrum Q ky 
is in reasonable agreement. One finds £ ml = 1.67 £ DNS 
and Q M4 = 1.14 Q DNS , whereas the values exhibit a 
substantial disagreement without a model, according to 
£ M " = 3.00 £ DNS and Q M0 = 1.42 Q DNS . 

In summary, the LES model leads to a far better agree- 
ment with the reference DNS than the runs without a 
model. As far as the overall heat flux levels (which are of 
prime importance) are concerned, the relative error with 
respect to the reference DNS is acceptable, amounting to 
less than 30% in all three cases considered. The model 
amplitudes c x and c y computed dynamically are found 
to be quite robust when varying the temperature gra- 
dient. The mean values are c x = 0.0155, c y = 0.0179 



in the weakly driven case, c x 
the CBC, and c x = 0.0140, c y 
driven case. 



0.0140, c y = 0.0212 for 
= 0.0219 for the strongly 



C. Robustness while varying the magnetic shear 

Next, we would like to investigate the robustness 
of the LES approach with respect to variations of 
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FIG. 5. Wavenumber spectra £ ky (at top) and Q ky (at bot- 
tom): Comparison between DNS and LES for the case of 
weakly driven ITG turbulence at am = 6.0. 



the magnetic shear s. The effects of the latter 
on plasma microturbulcncc has been the subject of 
numerous experimental^—, theoretical^, as well as 
numerical^r— studies. In this context, it was also 
found that negative magnetic shear can help improve the 
plasma confinement in a tokamak by decreasing the level 
of turbulence. Apart from the CBC case, three highly 
resolved DNS runs have been performed, corresponding 
to reversed (s = —0.4), low (s = 0.2), or high (s = 1.4) 
magnetic shear cases, compared to the CBC standard 
value (s = 0.796). The DNS perpendicular grid is kept 
fixed compared to previous sections: 
while other parameters are those of the CBC. 

In the case of reversed shear, the free energy and free 
energy injection spectra peak at a slightly higher k y value 
(kyPi ~ 0.3 compared to k v pt ~ 0.2 for CBC), as shown 
in Fig. [7] The total free energy is very small compared to 
the CBC, indicating a low level of turbulence. This effect 
has already been observed in a previous numerical study 
based on the spectral heat flu?*^. The LES offers a satis- 
fying agreement with the reference DNS spectra, except 



N x x N y = 128 x 64, 
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FIG. 6. Wavenumber spectra £ ky (at top) and Q ky (at bot- FIG. 7. Wavenumber spectra £ ky (at top) and Q ky (at bot- 
tom): Comparison between DNS and LES for the case of torn): Comparison between DNS and LES for the case of 
strongly driven ITG turbulence at wt. = 8.0. reversed shear (s = —0.4). 



for an underestimation of the free energy injection peak. 
The total free energy agrees reasonably well with the ref- 
erence value, £ Mi = 1.35 £ DNS , while £ M0 = 2.12 £ DNS . 
Considering the heat fluxes, one finds the same trend: 
Q MJ = 1.04 Q DNS and Q M0 = 1.31 Q DNS . 

In Fig. [8]), the results for the case of low shear (s = 0.2) 
are presented. Here, the turbulence level lies between 
those of the reversed shear and CBC cases. The free 
energy and free energy injection spectra, £ ky and Q ky , 
are fairly extended, up to k y pi ~ 0.4. The use of a 
LES model prevents the accumulation of free energy at 
small scales, while it moderates the appearance of large- 
scale structures without suppressing them completely (at 
small non-zero k y ). The total free energy obtained by the 
LES model is a bit far from the reference DNS value, 
£M4 = i.78£dn Sj but much better than the estimate 
obtained without a model, £ M0 = 2.89 £ DNS . The dis- 
agreement regarding the heat fluxes is again found to 
be more acceptable, according to Q M4 = 1.17 Q DNS and 
Q M0 = 1.56 Q DNS . 

Finally, the results for the case of high shear (s = 1.4) 



are displayed in Fig. |H1 Although the turbulence level 
is slightly lower than for CBC parameters, the free en- 
ergy and free energy injection spectra very similar to 
the CBC ones. A very satisfying agreement between 
LES and DNS is found regarding both the total free en- 
ergy (£ M4 = 1.23 £ DNS ) and the total heat flux (Q M4 = 
1.04 Q DNS ). Again, without using a model the accumula- 
tion of free energy at small scales leads to larger differ- 
ences: £ M0 = 1.98 £ DNS and Q M0 = 1.15 Q DNS . 

In all three cases, an important consequence of the use 
of a LES model is that it prevents the accumulation of 
free energy at small scales. In addition, for low magnetic 
shear, unphysical features at small k y (which may de- 
velop due to the relative coarseness of the chosen grid) 
are strongly reduced. Although not shown explicitly, the 
k x spectra of free energy £ k * and free energy injection 
Q k * are always in good agreement. One notes that the 
total free energy appears to be a very sensitive diagnos- 
tic. The presence of a model considerably enhance the 
agreement between DNS and coarser simulations. The 
total heat flux is estimated with an encouraging relative 
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FIG. 8. Wavenumber spectra £ ky (at top) and Q ky (at bot- 
tom): Comparison between DNS and LES for the case of low 
shear (s — 0.2). 



FIG. 9. Wavenumber spectra £ ky (at top) and Q ky (at bot- 
tom): Comparison between DNS and LES for the case of high 
shear (s = 1.4). 



error of less than 20% for all cases. 



In contrast to the findings of the temperature gradient 
scan, the LES model amplitudes exhibit substantial vari- 
ations for changes in the magnetic shear. In particular, 
the value of c x is found to be close to zero in the reversed 
shear case (c x k, 5 • 10 -5 ) as well as in the low shear case 
(c x « 2 ■ 10 -3 ). Meanwhile, it departs from the CBC 
value (c x = 0.0140) only moderately in the high shear 
case (c x = 0.0102). The time traces of c x are shown in 
Fig. [TO] On the other hand, the c y values do not vary 
much; one obtains c y = 0.0192, c y = 0.0223, c y = 0.0212, 
and c y = 0.0185 for the four values of magnetic shear (in 
increasing order). These results are a reflection of the 
effect of the magnetic shear on the turbulence, which in- 
cludes the twisting of the perpendicular eddies along the 
magnetic field. 
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FIG. 10. Model amplitude c x as a function of time for different 
values of the magnetic shear s. 
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V. DISCUSSION 

In the present paper, a dynamic LES procedure has 
been applied to gyrokinetic turbulence as described by 
the GENE code. This approach provides an automatic 
calibration of the free parameters associated with dis- 
sipative GyroLES models. The dynamic procedure has 
been found to be robust in a wide parameter range of 
the logarithmic temperature gradient and the magnetic 
shear. 

Comparisons between DNS and GyroLES simulations 
have been based on free energy and free energy injection 
spectra, £ ky and Q ky . Generally, the use of a LES model 
has prevented the accumulation of free energy at small 
scales. While the differences regarding the total free en- 
ergy can exceed 50%, simulations without a LES model 
are even much more inaccurate, exhibiting relative errors 
up to about 200%. Moreover, when considering the total 
heat fluxes, the GyroLES results are really encouraging, 
with relative errors below about 20%. In terms of com- 
putational cost, the GyroLES approach has been found 
to save a factor of about 20, requiring only about 250 
CPU-hours per single run. This allows for the possibil- 
ity to perform nonlinear gyrokinetic simulations an any 
modest cluster with relatively little effort. 

Obviously, future GyroLES studies will have to take 
kinetic electrons and their contribution to the overall en- 
ergetics into account. It may be expected that GyroLES 
will also be of great benefit in this wider context, again 
leading to major savings of computer resources. Thus, 
GyroLES is likely to enable large parameters scans of 
gyrokinetic turbulence which can be used, e.g., to effi- 
ciently couple turbulence and transport codes (see, e.g., 
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